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Abstract 

We establish an extended version of the Einstein - Maxwell - axion model by introducing 
into the Lagrangian cross-terms, which contain the gradient four-vector of the pseudoscalar 
(axion) field in convolution with the Maxwell tensor. The gradient model of the axion- 
photon coupling is applied to cosmology: we analyze the Bianchi-I type Universe with an 
initial magnetic field, electric field induced by the axion-photon interaction, cosmological 
constant and dark matter, which is described in terms of the pseudoscalar (axion) field. 
Analytical, qualitative and numerical results are presented in detail for two distinguished 
epochs: first, for the early Universe with magnetic field domination; second, for the stage of 
late-time accelerated expansion. 
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1 Introduction 



Basic elements of the modern theory of the axion - photon (pseudoscalar-photon) coupling 
have been formulated thirty four years ago, in 1977-1978. Peccei and Quinn pQ in the 
context of the strong CP- problem have opened the discussion about new pseudoscalar 
particles, which could appear as a result of a spontaneous phase transition (PQ- symmetry 
breakdown in modern terminology). Due to the works of Weinberg [2] and Wilczek [3] now 
we associate this problem with the theory of pseudo-Goldstone bosons and with the so- 
called WIMPs (Weakly Interacting Massive Particles). For particles of this type the term 
axions was introduced by Wilczek; one can find the history of axions, e.g., in [I]-[BJ. At 
the same time, in 1977, studying the Equivalence Principle, Ni j9] suggested to describe 
the pseudoscalar - photon interactions on the electrodynamic language: the cross-invariant 
(j)F£ in F mn was introduced into the Lagrangian, where 4>, the dimensionless pseudoscalar field, 
appears in the product with the convolution of Maxwell tensor F mn and its dual F^ n . Later, 
on the basis of these ideas, new theoretical models were established and developed, namely, 
axion electrodynamics (see, e.g., [10] -[15] ) and its covariant generalizations: the Einstein - 
Maxwell - (dilaton) - axion models (see, e.g., [16] -[23] for review and references). Nowadays, 
the axion electrodynamics provides a theoretical base for the experimental search for the 
axion-photon coupling (see, e.g., [2J] - [28] )• The most intriguing application of the axion 
theory is connected with the concept of dark matter, whose contribution into the Universe 
energy balance is estimated to be about 23 % (see, e.g., [5j |8J [291 EQl EH [32] . If the hypothesis, 
that the dark matter is composed of axions, will be proved, and about 23% of the Universe 
energy will be prescribed to the contribution from the so-called background (relic) axions, 
we could conclude that the axion-photon coupling is the important element in the hierarchy 
of cosmic interactions. 

By analogy with classical Maxwell- Faraday electrodynamics, the axion electrodynamics 
admits generalizations of different types. There are various reasons for the extension of this 
theory; for instance, the nonminimal extension of the Einstein - Maxwell - axion theory 
(see, [331 El]) was motivated, in particular, by the problems of correct estimation of the 
polarization rotation effect for the electromagnetic waves running through the regions with 
nonstationary and/or nonuniform gravitational field. Other arguments for the extension of 
the axion electrodynamics come from the analogy with electromagnetic theory of moving 
anisotropic inhomogeneous (polarizable, magnetizable, rheological) media [351 ESI E7]: in 
these models the tensor of linear response depends on the macroscopic velocity of the medium 
as a whole, and its derivatives (see, e.g., [381 El]), thus providing the axion system to possess 
the birefringence, dynamo-optical, etc. effects. 

The gradient extension of the Einstein - Maxwell - axion model presented in this work has 
the following argumentation. First of all, when the gradient four-vector of the pseudoscalar 
field, Vj0, is nonvanishing, the axion - photon system can be considered as an anisotropic 
one - axis quasimedium, with a director dti, normalized four- vector, which is proportional to 
Vj0. Anisotropy of continua is expected to be inherited by permittivity tensors [35], thus 
providing the birefringence effect. However, the standard axion electrodynamics with the 
cross-invariant <pF^ in F mn in the Lagrangian does not take into account the anisotropy of the 
constitutive equations, and this weak point of theory can be amended by the gradient type 
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extension. The second motive to insert the gradient four-vector Vj0 into the Lagrangian 
comes from the analogy with the so-called derivative coupling in the nonminimal theory of 
the scalar fields [40j EH [42] and of the Higgs multiplets j43j EH EH] . Let us remind that in 
[HUSHES] i n addition to the standard invariant £i?$ 2 with Ricci scalar R and scalar field $, 
the new term aR Vj3>Vfc<l>, with the Ricci tensor R lk , was inserted into the Lagrangian. In 
our model in addition to the term (j)F^ in F mn we introduce admissible cross-invariants, which 
contain the convolution F* k V k <f). The third argument to extend the theory of axion-photon 
coupling is connected with the attempts to detect experimentally the birefringence in the 
vacuum magnetic field (see, e.g., [21] - [28]): the standard vacuum axion electrodynamics 
does not admit birefringence, while its gradient-type modifications permit it. 

The paper is organized as follows. In Sec. [2] based on the Lagrange formalism we ob- 
tain self-consistent set of master equations for the electromagnetic (Sec. 12.21) . pseudoscalar 
(Sec. 12. 3p . gravitational fields (Sec. 12.4)) . and discuss general properties of the model. In 
Sec. 12.51 we discuss the problem of bilingual description of the axion-photon interactions, 
and constraint the choice of the phenomeno logical parameters of the gradient model using 
two exact solutions for the induced electric field. In Sec. [3]we consider obtained master equa- 
tions by the example of anisotropic cosmological model with initial magnetic field, axionically 
induced electric field, dark energy in A-term representation and dark matter described in 
terms of pseudoscalar (axion) field. The reduced system of master equations (Sec. 13. ip is 
presented in terms of associated dynamic system (Sec. 13.31) . The epochs of dark matter dom- 
ination (Sec. I3.4f) and of magnetic field domination (Sec. 13.5ft are distinguished and studied 
analytically and qualitatively; the model as a whole is analyzed numerically. 



2 Two- parameter gradient-type model 
of the axion-photon coupling 

2.1 General formalism 

The gradient-type extension of the axionic Einstein-Maxwell model is based on the La- 
grangian formalism, the action functional under consideration is of the form 



/( H I ) A 1 1 
d 4 xV^ [—^+-F mn F mn +-^ n F mn +^l [-V m 0V m 0+n<A 

+ h 1 F mn F mn V p ^^+h 2 F mp F m ^V g ^ . 



+ 



Here g is the determinant of the metric tensor go-, V m is a covariant derivative, R is the 
Ricci scalar, k=^- is the Einstein coupling constant, A is the cosmological constant. The 
Maxwell tensor F mn is given by 

Fmn = ^m-^-n V ' n-^m = 2Vj m ^4. n ] , (2) 

where A m is an electromagnetic potential four- vector; F* mn =^e mnpq F pq is the tensor dual 
to F pq ; e mnpg = -±=E mnpq is the Levi-Civita tensor, E mnpq is the absolutely antisymmetric 
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Levi-Civita symbol with £' 0123 =1. The dual Maxwell tensor satisfies the condition 

\/ k F* ik = . (3) 

The first term — corresponds to the standard Hilbert - Einstein Lagrangian; the second term 
is the standard Lagrangian for an electromagnetic field; the third term is the pseudoscalar 
- photon interaction Lagrangian; the fourth and fifth terms constitute the pseudoscalar 
Lagrangian. The symbol stands for a pseudoscalar field, this quantity being dimensionless. 
The axion field itself, $, is considered to be proportional to this quantity $ = ^/ (p with a 
constant \l/o- The function V(4> 2 ) describes the potential of the pseudoscalar field. Two last 
terms in the Lagrangian introduce new terms containing the gradients of the pseudoscalar 
field in the product and/or in convolution with the Maxwell tensor; the parameters Ai and 
A 2 are phenomenological coupling constants. 

2.2 Axion electrodynamics 

Variation of the action functional ([1]) with respect to the four-vector potential Ai gives the 
equations of axion electrodynamics 

V k H ik = , (4) 
where the excitation tensor H lk is given by the term 

H ik =F ik +<pF* ik +\ l F ik V q (f)V q <p+\2V [k <pF l]q V q (l). (5) 
Using the linear constitutive equations 

H ik = C ikmn F mn , (6) 
we readily obtain that the linear response tensor C tkmn takes now the form 

C ikmn = 1 ^ g im g kn_ g in g kmj V p </>V P ' 



It is well known that using the medium velocity four-vector U\ normalized by U i Ui = 1, one 
can decompose the linear response tensor C lkmn uniquely as 

^jikmn ^i\mjjn\jjk _|_ ^k\njj"m]jj-i | 

- \rf kl {^r 1 ) ls rf nns + rf kl U [m v l n] + rf mn U {i v { k] . (8) 

Here e im is the dielectric permittivity tensor, is the magnetic impermeability tensor, 

and v p m is the tensor of magneto-electric coefficients. These quantities are defined as follows 

Am c\/^iikmnjj tt 

e = zo UkU n 

v v m = VmkC tkmn U n = U k C mkln Vlnp . (9) 
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The tensors r] mn i and rf are the anti-symmetric ones orthogonal to U l and defined as 

Vmni = e mnls U s , V ikl = e Ms U s . (10) 
They are connected by the useful identity 

- V^Vmnp = €nsUlU S = A^ - A^A* , (11) 

where the symmetric projection tensor A lk is 

A ik = g ik _ jj^k _ ^ 

The tensors £«. and (/x -1 )^ are symmetric, but vik is in general non- symmetric. Clearly, 
these three tensors are orthogonal to U l . Using the expression (J7J) one can calculate tensors 
e im , (/i _1 )i m and v vra explicitly, yielding 



A im {V<t)f +V>V r 



U- 1 ). = A im [l + A 1 V 9 0V 9 0] + ^A 2 



(13) 



(14) 



(15) 



The symmetric tensors e lTn and (/i -1 )i m contain both new coupling parameters Ai and A 2 , 
while the cross-tensor vik, which describes optical activity effects (see, e.g., |16]), contains 
the second coupling parameter only. Here we used the standard definitions 



V = U m V r 



V fc = ATV r 



(16) 



for the convective derivative and pure spatial gradient, respectively. 

One can mention two interesting details. First, the tensors e im and (yU _1 )j m become 

spatially anisotropic, when the pseudoscalar field is inhomogeneous and A 2 7^ 0; this means 

that the axion-photon interaction of the gradient type can produce the phenomenon of 

birefringence in the course of electromagnetic wave propagation. Second, the cross-tensor 

u pm con t a i ns both symmetric and skew-symmetric terms (see the terms with A pm and rj pmh , 

respectively). It is well-known that the non- vanishing cross-tensor v %k indicates that the 

axionic vacuum as a (quasi) medium is optically active, and the rotation of the Faraday type 

takes place in the course of electromagnetic wave propagation. Thus, one can see directly 

from (j!5p . that the interaction of a new gradient type between electromagnetic and axion 

fields produces optical activity of new type, if the pseudoscalar field is non-stationary, i.e., 

± 

V(p 7^ 0, and inhomogeneous, i.e., 7^ 0, simultaneously. 
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2.3 Pseudoscalar field evolution 



Variation of the action functional (CQ) with respect to the pseudoscalar field (ft gives the 
equation 

V 9 [($«-e«) V p 0] +<PV> (0 2 ) = -J-F^F*™ , (17) 



where the tensor Q pq is defined as 



1 



® m = ^ [A i g pg F mn F mn +X 2 F mp F m 9 ] , (18) 

and the prime denotes the derivative with respect to the argument 2 . The pseudoscalar 
F^ nn F mn plays the role of the local source of the axion field, as in the standard theory of the 
axion-photon interaction. The novelty of the model is connected with the tensor 

g pq = (g pq - ® pq ) , (19) 

which in fact plays the role of an effective metric for the pseudoscalar waves: these waves 
propagate in axionic (quasi) vacuum similarly to the waves in the medium, so that the 
velocity of pseudoscalar waves depends on the structure of the electromagnetic field (see, 
e.g., [SI HU US] , in which color-acoustic analogs of these waves have been studied). 



2.4 Gravity field equations 

Modified Einstein's equations obtained by the variation of the action functional (pQ) with 
respect to the metric g pq can be presented in the form 

R pq - l -Rg pq = Ag pg + k fzjf> + T#> + \ X T$ + \ 2 T^f 

The stress-energy tensor of the electromagnetic field 

1 



pq 



rp(EM) 

pq 



7 gpq-FmnF 



F F ' 

1 ?m J q 



and the stress-energy tensor of the pure axionic field 



are presented by the well-known terms. Two terms 



Tff = --F mn F m "V p 0V g + T p (f M )V n 0V^ 



and 



T (2) 



pq 



F ql F v m V m (l>+F ml F™V p( l)+F ml F™V^ 



(20) 



(21) 



(22) 



(23) 



(24) 



describe new source-terms in the right-hand side of the gravity field equations. 

Thus, the complete self-consistent system of master equations describing the extended 
Einstein - Maxwell - axion model is derived: it contains the electrodynamic equations (see 
Q, ©), ©, the equation for the pseudoscalar (axion) field (see (fT7j) . ffTH]) ) and the gravity 
field equations (jUD-flMp. 
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2.5 Dark matter and alternative description 
of the axionic source-term 

Axions are considered to be candidates for the cosmic substrate indicated as dark matter 
(see, e.g., |S1 [2QJ 1301 E] ) • There are two different approaches to the description of axions: the 
first one is based on the pseudoscalar field theory, the second approach uses the terminology 
of the theory of continuous media. In the first approach the stress-energy tensor of the 
axions is quadratic in the gradient four-vector and is the direct result of the variation 
procedure. In the second case we deal with some phenomenological representation, and the 
best known construction is the stress-energy tensor of an imperfect fluid with the energy 
density W, heat flux four-vector q l and the (anisotropic) pressure tensor Vik- Let us assume 
that these quantities coincide, i.e., 

1 



% V^V fe 0-- to V m 0V m 0-^(r) 



WUiU k + qi U k + q k Ui + V lk = 4 A) . (25) 

Clearly, the compatibility of such equality is connected with some physical and mathematical 
constraints. Let us start with the discussion of the heat flux four-vector. Usually, the axionic 
systems, considered as a fluid, are characterized by vanishing heat flux four-vector 



q L = U p T^A iq = . 



(26) 



Mathematically, this can be obtained when the macroscopic velocity four-vector U l is the 

eigen- vector of the stress-energy tensor (the Landau-Lifshitz definition). Direct calculation 

± 

of q l using the left-hand side of the equality (T25]) yields q t =^QD<f)V t <j)=0. When the pseu- 

± . 

doscalar field is homogeneous, i.e., V l 0=O, or when the pseudoscalar field is stationary, i.e., 
V<ft=0, one obtains readily that q l —0. To provide the equality q l —0 in general case (inhomo- 
geneous and non-stationary simultaneously) one can assume that the macroscopic velocity 
four-vector is chosen to be orthogonal to the gradient four - vector of the pseudoscalar 
field, i.e., U<p=U l Vi<j)=0. This equation can be generally resolved, when Vk4> is a spacelike 
four-vector. The energy- density scalar W can be written as follows 



W = U p T^U q = -^l 



-pq 



(V(j)) 2 +V- V>V; 



(27) 



This quantity is always positively defined. The structure of the pressure tensor also attracts 
special attention. In many works the dark matter is considered to have the isotropic (Pascal) 
pressure P, moreover, in the cold dark matter model this quantity is vanishing P — > 0. The 
effective pressure tensor presented in terms of derivatives of the pseudoscalar field has in 
general case both isotropic and anisotropic parts 



ik 



1 
2 



V-(P0) 2 -V m 0V r 



(2* 
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Clearly, when the pseudoscalar field is homogeneous, the effective pressure tensor is spatially 
isotropic. Thus, keeping in mind the bilingual description of dark matter axions in terms of 
perfect fluid and in terms of pseudoscalar field, we have to provide two constraints on the 
pseudoscalar field model: first, the heat flux to be vanishing, second, the pressure tensor to 
be isotropic. Below we remind two explicit examples as illustrations that such models exist. 



2.5.1 First explicit example: spatially homogenous cosmological models 



Anisotropic spatially homogeneous cosmological model of the Bianchi-I type is known to be 
characterized by the metric 



ds 2 = dt 2 -a 2 {t) (dx 1 ) 2 -^) (dx 2 ) 2 -c 2 (t) (dx 3 ) 2 . 



(29) 



In the framework of spatially homogeneous Bianchi-I cosmological models, as well as for the 
isotropic FLRW models obtained by a=b=c, one assumes that all functions of state depend 

on time only, so that Vj0=O. The macroscopic velocity four-vector is of the form U 1 =Sq, 
and thus T>cf) = <p, where <fi is the time derivative of the pseudoscalar (axion) field. In these 
models the heat flux four- vector vanishes. We use here the well-known procedure (see, e.g., 
for details): the energy density of the dark matter can be expressed as follows 



<P 2 + V((P 2 ) =W{t). 



Clearly, the pressure tensor (|28|) is now isotropic 

where P is the Pascal (isotropic) pressure of the dark matter 

1 



2 



P-V(cj> 2 ) =P(t). 



These two relations yield the time derivative 



1 



i = ± _^( f)+ p( f ), 



and the pseudoscalar field itself 

- 0(t o ) = ±-L fdt'^wif) 

M/n Jtn v 



P(t>) 



(30) 



(31) 



(32) 



(33) 



(34) 



in terms of dark matter energy density and pressure. The pseudoscalar field potential V(4> 2 ) 
can be reconstructed using the relation 



^0 



(35) 



and the function (j)(t) obtained in 
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Let us mention that in present epoch of the Universe accelerated expansion the ground 
state of the dark matter is considered to be non-relativistic, i.e., P « W. More precisely, 
the CMD model (Cold Dark Matter) assumes that P=0 and one deals with a dust. In this 
case using ( )32l) one can conclude that ^(f> 2 =^V((f) 2 ) in line with the well-known classical virial 
theorem. Thus, in order to obtain the CDM model one should put P=0 into (l33l)- (|35l) . and 
replace the total energy density W(t) by the rest energy density pc 2 , where p{t) is the mass 
density of the dark matter dust. 

The dielectric permittivity tensor can be now reduced to the scalar function e(t) 

e m =A m e(t) , e(t) = 1 + (Ai + \\^\ <P 2 . (36) 



Similarly, the magnetic impermeability tensor has the form 

1 1 

W)' W) 



The cross-tensor is now pure symmetric 

v vm = -0A pm . (38) 

The tensors e im and (p~ l ) im are spatially isotropic, thus, the birefringence effect is absent, 
if we restrict ourselves by spatially homogeneous cosmological models. The refraction index 
n{t) can be introduced by the relation 

n 2 (t) = e(t)p(t) = ^ 2^_rv_ 3g 

W WPW l+A^ 2 V ; 

When Ai=A2=0, one obtains that n 2 =l, the phase velocity of electromagnetic waves f p h=^y 
coincides with speed of light in vacuum. In general case the phase velocity depends on time 
through the function <ft 2 . The phase velocity can take infinite value, when the denominator 
of fl39|) vanishes; it is equal to zero, when the numerator of fl39|) vanishes. The function n 2 (t) 
can be (in principle) negative, and in this (unlighted) epoch the electromagnetic waves do not 
propagate in the Universe, do not scan its internal structure and can not bring information 
to observers. 

The Bianchi-I model with pure magnetic field is well-studied (see, e.g., [Ml H9] for re- 
view). The novelty of our model is that due to the axion-photon coupling the magnetic field 
B 3 (t)=—-^Fi2 7^ produces the electric field, so we consider the electromagnetic configu- 
ration with supplementary component E 3 =F 30 (t) ^ 0. When all the quantities depend on 
time only, the equations ([3]) are satisfied, if Fi 2 =const. The equations (jl]) with the excitation 
tensor (jSJ) give the following exact solution for the electric field 

F 30 (t) _ F 12 <p(t) + const (4Q) 



a(t)b(t)c(t) [l+ (Ai+|A 2 ) <P 2 {t) 



Clearly, this electric field can reach infinite value, if the parameter Xi + ^X 2 is negative; 
this means that at some critical moment the invariants of the electromagnetic field can be 
singular. Thus, excluding the singularity of this type we require that the coupling constants 
satisfy the inequality Ai+|A 2 > 0. 
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2.5.2 Second explicit example: static models with spherical symmetry 



Static spherically symmetric configurations can be characterized by the metric 

ds 2 = a 2 Ndt 2 -^dr 2 - r 2 (d9 2 + sin 2 



(41) 



where the metric functions a(r) and N(r) depend on the radial variable r only. Assuming 
that also depends on r only we obtain readily that q % = 0, and 



W{r) 



2 

P± = vl 



V((f) 2 ) + N(f>' 2 
W{r) 



Pi = VI 



V{(t) 2 ) - Ncj)' 
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(42) 

(43) 
(44) 



The prime denotes the derivative with respect to r. The profile of the pseudoscalar field <p(r) 
can be now reconstructed using the relation 



1 



W(r) + Pu(r) 



(45) 



In order to find constraints for the coupling parameters Ai and A2, let us consider the model 
with magnetic monopole and radial electric field induced by axion-photon interactions. The 
Maxwell equations admit now the following solution 



dip 



\x sin 9 



?0r 



+ const 



(46) 



In order to guarantee that the solution is nonsingular for arbitrary pseudoscalar field we 
should assume that Ai+|A 2 < 0. Combining two constraints obtained for cosmological 
model and for this static model we see that only the requirement Ai+|A 2 =0 is appropriate. 
Below we consider the one-parameter models with \ 2 =—2\i. 



3 Cosmological application: 

the toy model of the Bianchi-I type 

The gradient-type extension of the Einstein - Maxwell - axion model formulated here attracts 
attention, since several new tests can be suggested on the base of this model. In this paper we 
consider shortly one cosmological application based on the anisotropic Bianchi-I type model. 
Our choice can be argued as follows. On the one hand, the model of anisotropic expanding 
Universe with a strong magnetic field and axionic dark matter is the best theoretical model 
for the study of the axion-photon coupling and its consequences. On the other hand, the 
impressive detailed analysis of the recent observations of the Cosmic Microwave Background 
(CMB) radiation (see, e.g., [50J) gives a unique possibility to fit free parameters of the model 
under consideration and to check the main consequences. Our strategy is the following: in 
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this paper we clarify the roles of guiding parameters of the model and illustrate (analytically, 
qualitatively and numerically) basic tendencies in the evolution of main state functions. Next 
paper, entirely devoted to the fitting of the guiding parameters of the model and based on 
the comparison with results of [50J, is now in preparation. 

Let us define concretely main problems, which we intend to consider in this cosmological 
application of the established model. First of all, we intend to demonstrate that the axion- 
photon coupling in presence of strong cosmological magnetic field can be a source of a 
substantial growth of the pseudoscalar (axion) field (f>; keeping in mind the hypothesis that 
the background (relic) axions form the dark matter, we hope to explain large contribution 
from this cosmic substrate into the total energy- density (23% in our late-time epoch). Second, 
since the axions provide spatially isotropic source of the gravitational field, we intend to show 
that the substantial growth of the number of axions, caused by the axion-photon coupling, 
could accelerate the isotropization process of the Universe. Third, we hope to illustrate the 
idea that the axion-photon coupling can be a reason for multi-stage inflationary-type regimes 
of evolution of the early Universe. 

Anisotropic cosmological Bianchi-I model characterized by the metric ( 129]) is assumed to 
be constructed using three constituents: the dark energy (it is presented by the A-term), the 
dark matter (it is described by the axions in terms of pseudoscalar field) and the electromag- 
netic field (it has originally the magnetic component, but the electric component appears 
due to the axion-photon coupling). The total set of master equations can be reduced now 
as follows. 



3.1 Reduced master equations 

3.1.1 Exact solution to the electrodynamic equations 

Due to the axion-photon coupling the magnetic field B 3 (t)=— ^-^12 7^ (Fi 2 =const) pro- 
duces the electric field E 3 =F 30 (t) 7^ 0. We consider the model with Xi+^X 2 =0, thus the 
equations (jl]) with the excitation tensor ([5]) give the following exact solution for the electric 
field 

F 30 (t)= /^ww (47) 
v ; a(t)b(t)c(t) v ; 

The corresponding physical components of the magnetic B(t) and electric E(t) fields, given 
by 

ms fi&5 = -^, (48) 



E(t)sy f=^ = ^. = Bmt ), (49) 

are proportional to the constant F 12 , which we will consider as a guiding parameter of the 
model. 
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3.1.2 Reduced equation for the pseudoscalar field 

The axion field function <p(t) can be found from the equation 



1 d 

abc dt 



abc<fi \ 



^l F 12 

a?b 2 



p2 N 

12 -^IV 



a 



2 b 2 -0' 



(50) 



When Ai is positive, there exists a moment t( cr it) , for which the coefficient in parentheses in 
the left-hand side of (!50|) takes zero value, and the equation becomes degenerated. In order 
to guarantee the absence of singularity in the equation (1501) . we assume below that Ai < 
and thus A 2 =— 2Ai > 0. 

3.1.3 Reduced Einstein equations 

Evolution of the gravitational field is described now by the following system of nonlinear 
equations: 



a b at be 

-T + + 7 - =A + K 

oo ac be 



F 2 



2^ V + ^ + 2aW 



(l - Ai0 2 + <p 2 ) 



(51) 



b c be 
b c be 



a c ac 
- + - + — 

a c ac 



A + K 

A + K 



a b ab 
-+-+--= A + K 
a b ab 



+ 



F 2 

r \2 

2a 2 b 2 
F 2 

r 12 

2a 2 b 2 



(l + Xi<p 2 + 



1 + Ai0 2 + <P' 



(52) 
(53) 
(54) 



Clearly, the sources in the right-hand sides of fl52|) and f|53|) coincide, since both the original 
magnetic and the induced electric fields are directed along one axis (Ox 3 ). This is a good 
motivation to use the model with the local rotation symmetry (LRS), which assumes that 
a(t)=b(t) c(t). Below we use this simplification and distinguish the longitudinal (along 
Ox 3 ) and transversal (on the plane x l 0x 2 ) dynamics; to describe them we introduce the 
longitudinal Hubble function Hn = -, and the transversal Hubble function H—-—\. 



3.2 Constant solution with hidden Aj and A2 

Let the potential of the axion field be quadratic in 0, i.e., V (<fi 2 )=m 2 <fi 2 . Then for the special 
set of the guiding parameters of the model the total system of equations (!50l) - (l54l) admits 
the exact solution, for which the pseudoscalar field is constant 4>(t)—(f>o 7^ 0, and the scale 
factors a and b are constant, i.e., a(t)=ao, b(t)=bo, the only factor c(t) being the function of 
time. Indeed, the equation (15T)1) is satisfied if 



FI2 



m 2 a\b 2 ^l . 



(55) 



12 



In this case the equation (!5T|) determines the value 0o by the relation 

This is possible only for the anti-de Sitter-type model, when the cosmological constant is 
negative, and 

|A| > ^K,m 2 y 2 Q . (57) 

Because of relation (!55l) the equation f|54|) turns into identity, and the equations (1521) and 
(JS3J) yield the key equation for the scale factor c(t) 

- = -n 2 , n 2 = |A| + -Km 2 ^l . (58) 
c 2 

Let us emphasize that the parameters Ai and A2 are hidden in this case: they do not appear 
in the equations since in all the key equations they are coefficients in front of the time 
derivative 6. The solution to ([58]) 



c(t) = c(0)cosfit + ^ sin fit, (59) 

describes harmonic oscillations of the Universe in the direction along the magnetic and 
electric fields. 

3.3 Associated dynamic system 

Let us consider the Bianchi-I model with LRS and introduce the following convenient di- 
mensionless variables 

x =4fL, - = xH-, (60) 
a(to) ' dt dx 1 

d ^ 7 H " 

Tx^ Z = H$ 

and auxiliary parameters 

iy 2 = -X 1 H 2 >0, H =H(to), (62) 

kF 2 

a = K*l, A = XH 2 , 0= I 2 . (63) 

n a {t ) 

For the potential of the pseudoscalar field we use here the (simplest) linear ansatz 

n^ 2 ) = tf 2 (7 + /iV) , (64) 

and formally put 7=0, since this constant can be included into an effective dimensionless 
cosmological constant A. In this terms the whole set of equations splits into one equation de- 
scribing the longitudinal dynamics and three-dimensional dynamic system for the description 
of the transversal dynamics. 



X = <f>(x), Y = x—J, Z = — 21 (61) 
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Remark about initial data 

The moment t=to of cosmological time, for which x=l, relates to some intermediate point, 
which splits the Universe history into two intervals: the second one, 1 < x < oo, finishes 
by the so-called late-time period, the first interval, < x < 1 starts with the so-called early 
Universe. Clearly, the replacement x — > 1/x converts the last interval into 1 < x < oo and 
effectively changes the time arrow. We choose the moment to so that at t=t Q the effective 
contribution of the magnetic field into the total energy of the Universe substratum is equal 
to the contribution of the dark matter. More precisely, we require that 

m 2 ^ 2 a\t ) = F 2 2 , (65) 

and thus (3=afi 2 (let us remind that vn}=^ 2 H'^). The parameter F 12 related to the initial 
value of the magnetic field is not fixed, it is a fitting parameter. In other words, the interval 

< x < 1 can be indicated as an epoch of magnetic field domination, while the interval 

1 < x < oo is in these terms the epoch of dark matter domination. The initial data for other 
unknown functions, namely, 4>{t ), <fi(t ), a(t ), a(t ), c(t ), c(t ), being necessary elements 
of the model, play different roles in the numerical and qualitative analysis. As we have 
mentioned above, the moment t=t relates to x—1, thus the values (fi(t ), <fi(t ) are codified, 
in fact, in the parameters X(l) and Y(l); we listed these parameters in the captions for 
the Fig.l, Fig.3, Fig.5. The values a(t Q ) and c(i ) are included into H(t )=a(t ) /a(t ) and 
H\\(to)=c(to)/c(to), respectively. The initial value of the transversal Hubble function H(t ) 
is not fixed, it is a fitting parameter. We analyzed numerically the ratio H(t)/H(to), so that 
the initial value of this function is equal to one (see the finishing point of the curve displayed 
on the panel (a) of the Fig.4, and the starting point of the curve on the panel (a) of the 
Fig.3). In our computational scheme the initial value of the ratio H\\(to) / H (to) (according to 
f l5~T|) ) is not an arbitrary parameter, it is predetermined by the initial values <p(to), <fi(to) and 
by the values of the guiding parameters a, z/, /i (see the right-hand side of floTj) ). Mentioned 
parameters are listed in the captions to the Fig.l, Fig.3 and Fig.5; the corresponding initial 
value for the ratio H\\(to) / H(to) can be found from the plots. Finally, the initial values of 
the scale factors, a(t ) and c(t ), are hidden parameters in our computational scheme, they 
also can be considered as fitting parameters. 



3.3.1 Evolution of the longitudinal scale factor 

The longitudinal Hubble function H\\=- can be found from the equation ( 15T|) (with a=b). 
In terms of variable x the equation for c(x) reads 



d . V c 
X Tx hgC= 2Z 



(66) 



where V c (x) is given by the function 



V C =X-Z+- 



// 



X 2 +ZY 2 +t-- (l+u 2 ZY 2 +X' 

x A v 



(67) 



The dimensionless function, which we are interested to find, is j^ = ^^ 
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3.3.2 Transversal 3D-dynamic system 

The dimensionless functions X, Y and Z, introduced by (!6lil . are coupled by the set of 
equations 

x-^-X = Y, X -^Y = V Y , x-^-Z = V z + X, (68) 

\JjiAj \AJtXj \AjJu 

where Vy and Vz are the following nonlinear functions of variables x, X, Y, Z, and of the 
guiding parameters u, /i, a, A: 



y Z(:e 4 +z/ 2 /i 2 ) 



Vc+Vz (x^V) 
2Z (x 4 +z/> 2 ) 

,2 



(69) 



y z =| (^X 2 -^ 2 ) -3^+^ (l-z/ 2 ZF 2 +X 2 ) . (70) 
Below we analyze this dynamic system qualitatively and numerically. 

3.3.3 Numerical analysis: 

general properties of magnetic, electric and pseudoscalar fields 

Numerical analysis of the model is fulfilled for various values of the dimensionless guiding 
parameters v,fi,a,\ and initial values X(1),Y(1). Below we discuss main tendencies and 
typical results using the most illustrative examples. Let us mention that the magnetic field 
B(x) = ^r- (see (jlSjl ) decreases monotonically, and its profile does not require illustrations. 
The electric field, induced by the axion - photon interactions, E(x) = B(x)(j)(x), (see (jUj) ), 
in fact inherits the properties of the axion field <f)(x), and we do not display the plot of its 
profile. Typical behavior of the pseudoscalar field <f>(x) (obtained from the solution for X(x)) 
is displayed on the Fig.l for the whole interval < x < oo. The line number 2 presents the 
plot 4>(x) for the model with Ai=A2=0, i.e., for the case when the gradient-type terms in the 
Lagrangian (TI)) are absent. Clearly, late-time evolution of the pseudoscalar field does not 
feel the presence or absence of the terms with Ai and A2; the interactions described by these 
terms are important in the early Universe. It is important to stress that when Ai=A2=0, 
the modulus \<p\ is infinite at x — > (see the second line), while for gradient-type models 
the pseudoscalar field typically tends to finite value at x — > 0. Another typical feature is 
a non-monotonic character of the 4>(x) profile: the pseudoscalar (axion) field grows in the 
early Universe, reaches the maximum value 0( max )(z/, fi, a, A) (at x > 1), then decreases at 
x — y 00 in a quasi-oscillatory regime. 

Concerning the plots for the longitudinal Hu—- and transversal H-- Hubble functions, 
below we consider them separately for the intervals < x < 1 and 1 < x < 00, since the 
appropriate visualization of the profile details requires to use different scales for these plots. 



3.4 The epoch of dark matter domination (x > 1): 

Asymptotic behavior and late-time isotropization 

The terms containing Fl 2 ar' 2 b~' 2 in the key equations of the model become vanishing at 
a — > 00, b(t) — > 00 (see, e.g., (I50l)-(l54l)). Both remaining sources of the gravity field, the 
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Figure 1: The plots of the pseudoscalar field, <j){x). The first line relates to a=0.1, z/ 2 =l, A=l, 
/z 2 =0.1, X(l)=0.5, y(l)=0.5. The example 2 differs from the example 1 by the condition 
Ai=0=A2, i.e., v 2 =Q. The line 3, with the guiding parameters a=0.01, z/ 2 =3, A=l, /i 2 =100, 
X(l)=0.5, y(l)=l, demonstrates a quasi-oscillatory regime of the asymptotic behavior at 
x — > oo. 

axionic dark matter and dark energy (A term) are spatially isotropic. Comparing the left- 
hand sides of the equations (i52)l - (1q"4"I) . and using the LRS ansatz, one can show that the 
appropriate solution in this limit is the model with In other words, the scale factors 

a=b and c become proportional to one another, thus, with rescaling of the coordinate x 3 , 
we can obtain the asymptotically isotropic model with a(t)=b(t)=c(t) — » oo. In this limit 
we deal with standard isotropic homogeneous FLRW-type model, in which the cosmological 
constant and the axionic dark matter predetermine the Universe evolution. 

3.4.1 Solutions with de Sitter asymptote 

The late-time isotropization at A > can be studied using the following qualitative argu- 
ments. At x — > oo the associated dynamic system becomes autonomous 



— X = Y, r = log x , 



(71) 




(72) 




(73) 




(74) 
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3]) has only one stationary point 

A 



The autonomous subsystem (1711 

X* = , K* = , Z. 
In the vicinity of this point the equation ( 1731) yields 



A 



Z(x) ^- + C 3 x- 
and we obtain the de Sitter asymptote 



I A 



C» x- 3 



— i , a(t) — >■ a(to) exp 



(75) 



(76) 



(77) 



The phase plane XOY describes (in the appropriate dimensionless variables) the behavior of 
4> and cf). It is easy to show that the stationary point (0, 0) on this phase plane is a stable 
node, when A > |m 2 , or a stable focus, when A < |m 2 . In both cases the pseudoscalar 
field itself and its time derivative vanish exponentially with time. Fig.2 displays the example 
of a phase portrait of this 2D dynamic system. 




Figure 2: Phase portrait of 2D dynamic system at x — > oo reconstructed for a— 2, A=l, 
(1=1. The stable stationary point corresponds to X=0 and Y=0, i.e., 0(oo)=O and </>(oo)=0. 

Concerning the longitudinal dynamics, in the vicinity of the stationary point the scale 
factor c(t) behaves as follows 



d_ 

dr 



log c — > 1 , c(t) —> c(t ) exp 



it -t ) 



(78) 



displaying explicitly the isotropization law. 
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3.4.2 Special solution with anti-de Sitter asymptote 

When A is negative the special solution of the isotropic model exists, for which the pseu- 
doscalar potential is linear in time. Using the simplest potential V((j) 2 )=m 2 <j) 2 we can reduce 
the master equations for the isotropic stage to the pair of equations containing, first, the 
standard equation 

4> + 3H<p + m 2 = (79) 
for the pseudoscalar (axion) field, second, the equation 

3H 2 = A + -K^ 2 (m 2 4> 2 + 2 ) (80) 

for the Hubble function H(t) = -. We would like here to attract the attention to the specific 
exact solutions linear in time 









, t* — 





0(t) = 0(0) ( 1 - f ) , U = k\ ^*o0(O) , (81) 



im 2 K 



ff(*) = y-g-* r o0(O)(l--J , (82) 

which hold, if A=— |m 2 . At the moment t = t* the pseudoscalar field and the Hubble 
function change their signs. This solution gives the Gaussian function for the scale factor 

a{t) = a{Q exp [~\rn 2 {t - t*) 2 J . (83) 

For the time interval t < t* this function grows, and the acceleration parameter 

. . d H 3 , 

-^— = 1 + - = !--^-^ (84) 

is positive, when t < t* — v3/m. Clearly, according to this model the accelerated expansion 
of the Universe will be changed by decelerated expansion, then the scale factor will reach 
the maximal value a( max ) = Q>(t*), and finally, the epoch of collapse will take place. 

3.4.3 Numerical analysis of the transversal and longitudinal Hubble functions 

at x > 1 

The results of numerical analysis of the dynamic system ( l66i) -( 170j) for the interval 1 < x < 
oo are presented in the Fig. 3. The panel (a) illustrates the behavior of the ratio H/Hq, 
obtained as square root of the function Z, where the transversal Hubble function is defined 
as H—-. Similarly, the panel (b) displays the behavior of the longitudinal Hubble function 
H\v=-. Both graphs illustrate the tendency to isotropization, i.e., both H and H\\ tend 



asymptotically to the value H oa =^j A/3. The type of graph is predetermined by the value 
of A: when A > 3 the plots start below the asymptote, when A < 3, they start above 
this asymptote; we used the second case for the illustration. Depending on the guiding 
parameters u, /i, a, A the profiles H(x) and Hu(x) can have zero, one, two, etc. extrema. The 
lines indicated by the number 3 relate to the model with Ai=0=A2- 
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Figure 3: Behavior of the Hubble functions in the epoch of dark matter domination (1 < 
x < oo). Panels (a) and (b) display the functions H/Hq and H\\/Hq, respectively. The lines 
indicated by the number 1 relate to a—1, v 2 =20, A=l, yU 2 =l, X(l)=0.5, F(l)=l; for the 
lines 2 the parameters a=0.05, z/ 2 =15, A=l, yU 2 =20, X(l)=0.5, Y(l) = l are used; the line 3 
illustrates the model with Ai=0=A2- 

3.5 The epoch of magnetic field domination (x < 1): 
Anisotropic early Universe 

When the values of the scale factors are small, the terms with F 2 2 dominate over the terms 
with \I/q in the key equations fl50|) - fl5^|) . We consider the dynamic system fl66|) - fl70|) numeri- 
cally in the interval < x < 1; the results are illustrated by the Fig. 4, Fig.5 and Fig.l. 

The distinguished feature of the graphs is the following: for various sets of guiding 
parameters there exists a pair of extremum points on the curve H(x)/H , the minimum (at 
x=X( m in)) an d the maximum (at x=X( max )). Changing the set of guiding parameters one can 
shift these extrema (see, e.g., first and second lines on the Fig.4). In the interval < x < 
£(min) the function H(x)/H grows infinitely at x — > 0. In the interval X( m in) < x < X( max ) 
the function H 2 (x) can be effectively fitted by the function H 2 oc x~ 4 , while in the interval 
i( miH ) < x < 1 it behaves as H 2 oc x~ 3 . The reconstruction of the vector field for the 
integral curves on the phase plane ZOY at small x (see Fig.5) confirms, that, indeed, the 
quantity Z=H 2 (x) / Hq decreases, then passes through the minimum and, finally, tends to 
infinity at x — > 0. Auxiliary dashed lines on the Fig.5 mark the positions of the minimum, 
clearly, the minimum shifts to smaller values of x when Y=xjj^ decreases. Such behavior 
can be explained qualitatively as follows. At x — y the dynamic system fl67|) - fl7D|) can be 
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Figure 4: Behavior of the Hubble functions in the epoch of magnetic field domination (0 < 
x < 1). The guiding parameters are the same as on the Fig.3. Additional window visualizes 
the behavior of the curve H(x)/H near the local minimum. 



transformed into 



d 

x- r X = Y, 

ax 



x—Y = Y 

ax 



4- 



2x 4 Z 



1+X' 



+ 



X 

U*z 



x—Z 

ax 



i i 

(XV u o 

3 + 7~Y 

2x 4 



+ 



2^ 



(l+X 2 ) 



and the approximate solutions have the form 



X ->x 



1- 



x 



2az/V(l+^o 2 ; 



, Y -> Y x b + 



2X n 



ai/V(l+*o) 



Z -> 



l+ia/i 2 (l+X 2 ) x- 3 -ia/i 2 (l+X 2 )x- 4 



(85) 
(86) 

(87) 



(89) 



Here X =X(0), and we assume that the fitted integral curve passes through the point 
Z(l)—1. The local maximum of this curve corresponds to the parameters 



^(max) 



2a/i 2 (l+X 2 ) 



i+W(i+x 2 ; 



; (max) " 



q/x 2 (l+X 2 ; 

6 X (max) 



(90) 



The formulas (j88p - (|9U|) give a nice qualitative illustration of the results obtained numerically. 
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Figure 5: Local fragment of the vector field picture for the integral curves on the phase 
plane ZOY at x — > 0. The point (0, 0) plays a role of a saddle point. Dashed lines indicate 
the points, which relate to the minimum positions of the curves depicted on the Fig. 4. 
The illustrative parameters are chosen as follows a— 15, z/ 2 =l, A=l, yU 2 =0.01, X(l)=0.5, 

y(i)=-2. 

4 Discussion 

Extended Einstein-Maxwell-axion model established in the work inherits main properties of 
the standard covariant model of pseudoscalar - photon coupling and introduces new interest- 
ing details; it is worth studying these new features in the cosmological context, taking into 
account (hypothetic) axionic nature of the dark matter. The analysis of the first subset of the 
total self-consistent system of master equations, namely, of the equations of extended axion 
electrodynamics (jlj),© an d (J3J), shows that in addition to the well-known magneto - electric 
effect the birefringence phenomenon is admissible. This phenomenon can be visible, when 
the gradient four-vector of the pseudoscalar field, has non-vanishing spatial part (see, 
e.g., the formulas ( TT3|) and (fl4"]l for the anisotropic permittivity tensors). In the framework 
of spatially homogeneous Bianchi-I cosmological model the birefringence does not appear, 
however, the magnetic field produces the electric field collinear to the magnetic one, due to 
the axion-photon interactions. This axionically induced electric field (see (j40p ) can be an in- 
teresting addition to the large-scale cosmological magnetic field, the cosmic substrate, which 
seems to be very important for the Universe structure and evolution (see, e.g., |5T| 152]). 

The extended equation of the pseudoscalar field evolution (see (IT7|) with (fTB]) ) differs 
from the standard equation by the propagation operator; this means, in particular, that the 
velocity of pseudoscalar waves propagation depends on the direction and intensity of electric 
and magnetic fields. In the cosmological context the most important consequence of the 
model extension for the axion field evolution is the following. The pseudoscalar source in the 
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right-hand side of (ITTj) being proportional to the scalar product of the electric and magnetic 
field vectors, is equal to zero, when the cosmological model includes the magnetic field only. 
Since due to the axion-photon interactions the initial magnetic field produces the electric 
field, this pseudoscalar source happens to be switched on in the early Universe and provides 
the (essential) growth of the axion field in the epoch of magnetic field domination (see, e.g., 
the Fig.l). One can expect that this (magneto-electric) mechanism of the pseudoscalar field 
production plays a nontrivial role in the process of the background (relic) axions creation, 
which is important for understanding of the dark matter problem. 

Extended gravity field equations (|20l)-(l24l) contain new source-terms, which modify the 
gravity field dynamics in the presence of the axion-photon interactions of the gradient type. 
In the cosmological context such modifications are shown to be essential, when the contri- 
bution from the magnetic field into the total cosmic stress-energy tensor dominates over the 
contribution of the dark matter, i.e., in the early Universe (see Fig.4). 

The gradient model of the axion-photon interaction is the two-parameter one, i.e., two 
coupling constants Ai and A2 have been phenomenologically introduced, when we formulated 
generic Lagrangian with gradient-type cross-invariants ([I]). Keeping in mind the principle, 
that it is preferable not to introduce "New Constants of Nature", we propose to use the 
following relationships. First of all, in order to guarantee that the electric field, produced by 
interaction of the initial magnetic field and axions, is nonsingular at least in two physically 
different situations (cosmological and static spherically symmetric models), we have shown 
that Ai and A2 should satisfy the linear relation Ai + |A2=0 (see ( l4T)i) and (1461)); thus the two- 
parameter model was reduced to the one-parameter. Similarly, the axion field is guaranteed 
to be nonsingular (see (I50|) ). when Ai is negative, i.e., Ai= — v 2 . There are few possibilities 
to express the parameter v 2 in terms of well-known physical constants. For instance, keeping 
in mind (16 9 p one can put v 2 =l/fi 2 . Since in this case the coupling constant is reciprocal 
to the axion mass in square, we have a direct analogy with the Drummond-Hathrell model 
of the nonminimal coupling of gravitational and electromagnetic fields [53] (let us remind 
that the mentioned coupling constant is reciprocal to the electron mass in square, or more 
precisely, to the square of the Compton radius). 

Another possibility to choose the constant v 2 by a physically motivated manner follows 
from the analysis of the Fig. 4(a): we can choose the parameter v 2 so that the moment 
£(min), at which the curve H(x) reaches the local minimum, will coincide with the Planck 
time t(pi an ck)- In this sense, the fragment of the curve H/Hq between the local minimum 
and maximum, displayed on the Fig.4(a), can be interpreted in terms of inflationary - type 
evolution of the anisotropic early Universe. According to ( 19 Op the maximum value of the 
transversal Hubble function iJ( max )=i? ^/ Z( max ) depends on the choice of the guiding pa- 
rameters au 2 , Xq and H . For the parameters listed in the caption the dashed curve has 
the maximum estimated by the value -ff( max ) — 10 3 -£^o; since -ff( max ) oc a^ ax ) (see (l90l ). this 
value can be made as large as necessary, by using small values of X( max ). It is important to 
emphasize, that the reference curve related to the model with Ai=A2=0 has no extrema and 
thus can not be used for modeling of the inflationary-type evolution. Clearly, in order to use 
the WMAP data [50] for the fitting of the parameters of our model, we should estimate the 
reference value -ff(i SO tro P ) at the moment t(i SO tro P ), which relates to the starting point of the 
isotropic expansion of the Universe. According to the Fig. 3, for the given guiding parameters 
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the isotropic epoch with H ~ H\\ starts at x ~ 10. We hope to present the results of the 
fitting procedure in the next work. 

Returning to the questions posed in the beginning of Section EJ one can conclude the 
following. First, indeed, the axion-photon coupling in the presence of a strong cosmological 
magnetic field can produce a substantial growth of the pseudoscalar (axion) field <fi (see Fig.l, 
the fragments of the curves for x < 1), thus providing a sensible number of cold axions in our 
(late-time) epoch. Second, the large number of axions produced in the epoch of the magnetic 
field domination can accelerate the process of the Universe isotropization (see Fig.3), and 
this isotropization process is multi-stage. 

We expect that the study of the gradient model of the axion-photon coupling will be 
interesting also in the context of astrophysical applications, and intend to consider static 
configurations of axionically coupled electric and magnetic fields in a future paper. 
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